# This file creates Figure SI.4

# install.packages("xtable")
# install.packages("ggplot2")
# install.packages("dplyr")

library(xtable)
library(ggplot2)
library(dplyr)

# load and rename mother estimates
het_res_mother <- 
  read.table("het_results_mother.txt")

# load and rename father estimates
het_res_father <- 
  read.table("het_results_father.txt")

plot_mother_fs <- 
  het_res_mother %>% 
  filter(V5 == 3) %>% 
  mutate(Mother = "Mother")

plot_father_fs <- 
  het_res_father %>% 
  filter(V5 == 3) %>% 
  mutate(Mother = "Father")

plot_fs <- 
  data.frame(rbind(cbind(plot_mother_fs),
                   cbind(plot_father_fs) ) )[, - 5]

colnames(plot_fs) <- 
  c("est", "se", "age", "year", "Mother")

## Create plots for first stage form

gg_fs <- 
  ggplot(plot_fs ,
         aes(y = est,
             x = age,
             ymin = est + qnorm(0.025) * se,
             ymax = est + qnorm(0.975) * se) ) +
  geom_errorbar(width = 0, size = 2) + 
  theme_classic() +
  facet_grid(Mother ~ year) + 
  geom_point(size = 4) + 
  scale_alpha_discrete(range  = c(0.5, 1),
                       labels = c("Father", "Mother"),
                       name   = "") + 
  ylab("Estimate") + 
  scale_y_continuous(limits = c(0,1.1),
                     breaks = 0:4/4) +
  xlab("Age of Firstborn (years)") +
  theme(axis.text = element_text(size = 14),
        axis.title = element_text(size = 16),
        legend.text = element_text(size = 12),
        legend.title = element_text(size = 12),
        title = element_text(size = 16),
        panel.spacing = unit(1, "lines")) 
